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Abstract. A granular assembly composed of a collection of identical grains may 
pack under different microscopic configurations with microscopic features that are 
sensitive to the preparation history. A given configuration may also change in 
response to external actions such as compression, shearing etc. We show using a 
mechanical response function method developed experimentally and numerically, that 
the macroscopic stress profiles are strongly dependent on these preparation procedures. 
These results were obtained for both two and three dimensions. The method reveals 
that, under a given preparation history, the macroscopic symmetries of the granular 
material is affected and in most cases significant departures from isotropy should be 
observed. This suggests a new path toward a non-intrusive test of granular material 
constitutive properties. 
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Particular attention has been devoted to the study of the static properties of 
granular assemblies over the past 15 years jl j. In that regard, a now famous and 
paradigmatic example is that of the stress distribution under a pile of sand. Depending 
on the manner this pile is built, different pressure profiles are measured. For example, 
when the material is poured from a point source, e.g. a hopper, then a stress 'dip' occurs 
at the base, under the apex of the pile 13 El , whereas when the material is poured 
in a rain-like procedure only a slight 'hump' occurs jl]. In this context, a discussion 
on the nature of the equations governing the distribution of stresses in granular media 
has emerged [3 13 El- The proposition of hyperbolic equations jH 13 UD] as opposed 
to the elliptic equations of elasticity lead to the need for newer more discriminating 
tests. These typically involved determinating the stress response to a localized force for 
a slab of grains (the response function). The initial motivation was that the hyperbolic 
equations predict a double peaked response (a ring in 3D) which would contrast with a 
single peaked profile expected for elliptic equations out of isotropic elasticity. It is now 
clear that for generic granular assemblies made of rough grains, an elastic description 
using elliptic equations is most likely correct [TTJ H2J EB] . But it is important to note 
that this need not be the case for samples prepared in the specific case of isostaticity 
(see [HI [13 H3 UZ1 HE] and references therein). 

One interesting outcome from the previous studies is the fact that, as the layer 
geometry is simpler than that of the pile, this stress response is a powerful tool to 
explore the effects of preparation procedures and details of the microstructure. The 
ultimate goal is to be able to link the response profiles to the internal texture of the 
packing. For instance, the microscopic texture of the packing is quite sensitive to the 
details of external actions such as pouring, shearing, etc... One would like to have a 
probe of this structure which determines the effective elastic properties. 

In this paper, we focus on the sensitivity of stress response profiles to the system 
preparation. First, we give a summary of recent experimental measurements on two- 
dimensional packings made of photoelastic grains [T2*] I19|. Then in the second section, we 
show series of experiments on three-dimensional piling fUEOj- In the third section we 
present new data obtained from numerical simulations of granular layers. All these stress 
profiles and the impact on future investigations are finally discussed in the conclusion. 

1. Two-dimensional photoelastic experiments 

Experiments were carried out on 2D granular systems consisting of disks and pentagonal 
particles that were cut from a photoelastic material. In these experiments, it was 
possible to measure forces on individual particles internally - hence the significant value 
of the approach. In principle, it is possible to directly obtain the forces at contacts by 
solving an inverse problem for the stress fields internal to the particles that satisfy the 
patterns seen in photoelastic images. This is a formidable task to carry out for large 
collections of particles. Hence, we chose a simpler approach in which we calibrated the 
mean force applied to a particle to the gradient of the image intensity associated with 
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Figure 1. Force distribution in a two-dimensional packing of photoelastic grains in 
response to a localized force at the top (the gravitational part has been subtracted). 
Darker zones indicate larger stresses. Chain-like structures are clearly visible. 




Figure 2. Stress response averaged over ~ 50 trials like that of figure ^ The shape 
of this response shows two lobes for a regular packing of circular monodisperse beads 
(left), but only one lobe (middle) when the layer is disordered (pentagonal particles). 
When the layer is sheared beforehand, the response is skewed in the direction of the 
shearing (right). The typical height of these pictures is 15 — 20 grain diameters. 

that particle. Specifically, we determine the mean force vs. G 2 , where G 2 is the integral 
of the square magnitude of the gradient of the image intensity, integrated over a particle. 
The justification for this approach is contained in [T^ H5]. 

We then carried out a series of experiments with the following general approach. 
We created nearly vertical stacks of particles that were supported very slightly by a 
vertical glass plate. With this arrangement, the friction force between particles and 
the plate was tiny. Necessarily, this approach caused some gravitational gradients in 
the sample, which we removed during the subtraction process described next. Typical 
sample sizes contained roughly 1000 particles and the sample was usually about twice 
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as wide as it was deep. Each experiment then consisted in the following: (1) a sample 
was prepared as above; (2) images with and without polarizer in place were obtained; 
(3) a point force was applied and a third image was obtained; (4) the point force was 
removed and a fourth image was obtained to see if there were shifts in the positions of 
the particles. Generally, we rejected an experiment if there were shifts in the particles 
positions between the first and fourth images. We then took the difference between the 
images with polarizers in places with and without the applied force. 

This then yielded an experimental approximation to the response function. It is not 
a priori manifest that there is a linear regime available in these experiments. However, 
we were able to identify such a regime in all cases, and we restricted our observations to 
this linear regime. An additional point is the fact that for any given set of macroscopic 
conditions, such as types of particles, friction, etc. any given realization was unique, as 
shown in figure ^ Only by obtaining an ensemble of responses was it possible to generate 
the mean response function. Throughout these experiments, we generated ensembles 
of typically 50 different realizations for each set of conditions. These data also yielded 
information on the variance of the distribution, and hence the range of local outcomes 
that one might expect for a given experiment. In these experiments, we considered 
a variety of different conditions. Specifically, we considered ordered hexagonal and 
square packing, for which we varied the particle friction coefficient and we considered 
completely disordered packing of pentagons. In the last case, we considered controlled 
amounts of disorder through bi-disperse packing of disks. In most of these experiments, 
we considered weakly compressed states (compressed due to gravity) and in one set of 
experiments, we explored the effect of modest amounts of shear. 

Here, we focus on only two issues, namely the role played by spatial order within 
the sample, and the effect of modest shear on force transmission. The effect of order is 
highlighted by parts (a) and (b) of figure EJ The left side of this figure shows the force 
response in a spatially order hexagonal packing of particles, and the middle panel shows 
data for pentagons. In the ordered case, transmission occurs chiefly along the lattice 
directions of the packing, with some friction-dependent widening with depth. This 
is consistent with the general predictions of anisotropic elasticity [T^l I2H 1221 12H] but 
experimental limitations could not rule out predictions of hyperbolic models for a weakly 
disordered system J21J. The case of the pentagons shows a single peak that broadens 
linearly with depth, consistent with an elastic-like picture. The linear broadening with 
depth is not consistent with a diffusive- like model, such as the g-model [213 EE] • The 
right panel of this figure shows the response to a normal force applied to a sample that 
has been subjected to a 5° simple shear deformation. The shear sets up a strong network 
of force chains that align along a direction of 45°. Indeed, along this direction, a force- 
force correlation function suggests long-range ordering, limited here by the sample size. 
And, interestingly, the force response is most strongly refocused to a direction of about 
22° with respect to the vertical direction. 

In these experiments, we exploited the power of birefringence in polymeric 
materials, to directly show the decisive influence of local granular organization on the 
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macroscopic mechanical properties of a granular piling. As evidenced by the mean 
response to a localized force, this method reveals a connection between local symmetries 
of the contact forces and macroscopic anisotropy. On the other hand, the birefringence 
method is restricted to 2D. Thus, it is important to see whether the effects displayed 
in 2D would persist for a more conventional granular packing such as sand. This is the 
purpose of the following section. 

2. Three-dimensional experiments 

A priori, the assessment of the stress response to a small localized force in 3D is a 
challenging experimental question. First, there is no easy way to directly probe stresses 
in the bulk. Second, the direct estimation of stresses at the boundaries, especially 
the bottom boundary, puts us in a situation where the measured stress signal would 
always be overwhelmed by the hydrostatic pressure. The method that we employ here 
to determine the small force response in the large hydrostatic background is based on 
a lock-in detection of the stress response to a slightly modulated signal added to the 
main localized force. This turns out to be extremely sensitive and able to discriminate 
the response part of the signal [TTJ EOJ- We also checked explicitly that the responses 
obtained are restricted to a linear regime. 

An interesting test of our method is provided by the study of a model elastic 
material such as gelatine with a known Poisson ratio v = 0.5. It is first important to 
note that the necessary presence of boundaries for all practical purpose is an element 
that needs to be considered explicitly in an elastic analysis. In the case of the stress 
profiles of an isotropic semi-infinite slab, the response function is solely dependent on 
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Figure 3. Stress response at the bottom of gelatine layers of various thicknesses. The 
two elastic profiles show the typical accessible range of shapes for the response of an 
isotropic layer of finite thickness on a perfectly rough bottom. The Poisson coefficient 
v = 0.27 gives the widest profile, while v = 0.5 gives the most narrow one 
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Figure 4. Response profiles for the vertical stress, a zz . The open circle data points 
were obtained on dense, compressed layers of grains (a), while the filled ones are from 
rather loose packings (b). We have plotted together data for layers whose thickness 
varies from h ~ 30 to h ~ 60 mm. We observe that the response function of a loose 
piling is much narrower than that of a compressed one. As for figure |3 the range of 
elastic predictions is shown by the two profiles v = 0.27 and v = 0.5. 

the spatial dimensions but this nice property is lost for finite depth elastic layers. For 
example, boundary roughness and Poisson ratio appear directly in the isotropic elasticity 
solution [20] • The gelatine results are plotted on figure El The data have been rescaled 
by the thickness, h, of the layer and by the amplitude of the extra force. The comparison 
to isotropic elasticity is excellent. Note, however, that the sensitivity of the response to 
the exact value of the Poisson coefficient is not particularly high; for instance, solutions 
with v = 0.27 and v = 0.5 show very similar shapes. 

We now turn to response functions obtained with different preparations of 
Fontainebleau sand (d ~ 300/im). To date, we have tested granular assemblies using the 
four following preparation methods: (i) a loose packing, (ii) a dense packing, both cases 
having the direction of gravity as a symmetry axis [2U| > (hi) a sheared granular assembly 
and, (iv) an horizontal slab prepared using successive avalanches j2Z|- We present in 
figures E] and |3] cases (i), (ii) and (iii). In all situations a linear elastic-like rescaling of 
the response with depth h was obtained - the scaling here refers to the fact that elastic 
stress fields have the form a a p = l/z p f(x/z) (p = 2 in 3D, p — 1 in 2D). 

Except for case (i), the predictions of isotropic elasticity are more or less inconsistent 
with our experimental results. For instance, the dense packing prepared by successive 
compression of deposited layers gives a response which is far too wide to be well fitted 
by such a model. For preparations obtained either in a shearing box or by avalanches 
- cases (iii) and (iv) - a shift of the response maximum with respect to the vertical 
direction below the applied force is present. It directly attests to the presence of an 
anisotropic medium with a symmetry axis other than the vertical direction. From all 
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Figure 5. Response profile for the vertical stress a zz at the bottom of a sand layer 
which has been uniformly sheared along the x axis. The different symbols are for layers 
of various thicknesses (50 to 100 mm) and shear strains (2° and 5°). The maximum of 
the profile is clearly shifted in the direction of shear. Isotropic elasticity is not able to 
reproduce such an asymmetrical shape of the profile with a vertical overload. 



these data, as in 2D, we conclude the existence of a macroscopic anisotropy induced 
by the pile preparation. Similar conclusion could be drawn for 3D response functions 
experiments [2EJ EH] or simulations and calculations performed on highly symmetrical 
crystalline packing O EUJ EH E2 E3 . 

Recently, there have been attempts to describe granular assemblies using 
orthotropic elasticity [21] where the direction of larger compression is taken as 
the stiff direction and thus the main symmetry axis of the model. This approach can 
provide a good qualitative account of all the experimental results but the number of 
parameters involved is too high to provide a decisive interpretation. Therefore, there is 
a crucial need to measure other quantities based on the response function method, such 
as shear stresses or responses to a non- vertical force as in [T^]. This is also the object of 
a numerical study in the following section. A large collection of these experimental data 
can provide severe tests of any proposed mechanical modelling of granular assemblies 
under quasistatic deformation. We pursue this issue further via a numerical study which 
we discuss in the following section. 

3. Numerical simulations 

In order to complement the experimental data presented above, we performed extensive 
2D simulations of various grain packings obtained with different preparation methods. 
The control of all the parameters of the simulations as well as the ability to measure 
both micro (grain size) and macro (system size) quantities ensures a useful feedback to 
the experiments. In addition, a numerical approach allows an easier exploration of a 
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larger range of preparation procedures. Similar numerical studies have been recently 
carried out, either on granular packings jHUEH] or on amorphous elastic bodies such as 
Lennard- Jones particles 

3.1. Simulation model and parameters 

The system consists of a set of iV polydisperse discs, with random radii homogeneously 
distributed between R m in and R ma x = 2Rmin- This system is simulated using a discrete 
element method - molecular dynamics (MD) with a third order predictor-corrector 
scheme |SZj. The rheology of contacts is that of Kelvin- Voigt. All grains interact via a 
linear elastic law and Coulomb friction when they are in contact. The normal contact 
force f n is related to the normal apparent interpenetration 5 of the contact as f n = k n - 5, 
where k n is a normal stiffness coefficient (5 > if a contact is present, 5 = if there is 
no contact). Note that this linear behaviour is consistent with Hertz law for discs. The 
tangential component ft of the contact force is proportional to the tangential elastic 
relative displacement, with a tangential stiffness coefficient k t . The Coulomb condition 
\ft\ < p-fn requires an incremental evaluation of ft every time step, which leads to some 
amount of slip each time one of the equalities f t = ±tif n is imposed. A normal viscous 
component opposing the relative normal motion of any pair of grains in contact is also 
added to the elastic force f n to obtain a critical damping of the dynamics. We chose a 
value /i = 0.5 for both the static and dynamic friction coefficients between the grains. 
All grains can be considered as quite rigid, k = k n /P = 1000, where k is a dimensionless 
parameter which expresses the level of contact deformation and where P is a pressure 
related to the weight of grains assembly. The ratio k t /k n was set to 0.5 or 0.7 [SB]- The 
bottom of the system is a horizontal line of similar but fixed grains, with intergranular 
distances chosen in order to avoid grains escaping from the system. We use horizontally 
periodic boundary conditions. Finally, we have carried out simulations on eight different 
samples of N = 3600 particles. 

We considered two different system preparations. The first is a rain-like (RL) 
procedure where all grains are initially placed at the nodes of a 1 to 4 aspect ratio grid. 
When the deposition starts, all the constraints are removed at the same time, letting 
the particles settle down under gravity with random initial velocities. In the second 
preparation scheme, each grain is deposited individually (grain by grain, or GG) with 
no initial velocity, by choosing a random initial position in contact with a grain already 
present at the outer surface. The next particle is deposited after 100 MD time steps if 
all grains already present have at least one contact. The equilibrium criteria consists 
of the five following tests which are applied after each period of 100 MD time steps: 
(1) the number of gained/lost contacts between particles during this period is zero, (2) 
similarly, the number of sliding contacts is also zero, (3) the integrated force measured 
at the bottom of the layer is equal to the sum of the weight of all the grains within a 
minimal tolerance, (4) all the particles have formed at least two contacts, and (5) the 
total kinetic energy is lower than some low threshold. 
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Figure 6. Contact and force angle distribution for the rain-like (left) and grain by 
grain deposition (right) procedures. Black lines are for the angles of contact whereas 
gray ones are for the orientation of the forces - all angles are measures with respect to 
the horizontal axis. The circles are guidelines indicating an isotropic distribution. The 
texture anisotropy of the grain-by-grain deposition is striking. These statistics have 
been obtained over about 50000 (RL) or 32000 (GG) contact pairs, and the bin width 
for each histogram corresponds to an angular interval of 5° 




Figure 7. Contact forces in response to an incremental small localized force at the 
top surface. Lines are bolder for larger forces. Black (red/grey) lines indicate contacts 
where the force has increased (decreased). The numerical results presented in this 
paper have been obtained with assemblies of 3600 particles forming typical layers 
of aspect ratio 1 x 4, of 21 grains thickness, and with horizontal periodic boundary 
conditions. 

We then carried out a microscopic characterization of the deposited layers and 
computed such quantities as the contact and force angle distributions. These are plotted 
in polar representation in figure El It is striking to see the difference between the two 
plots: the first one is close to a circle (quasi-isotropic distribution) whereas the second 
one shows four distinct lobes. Similar features have been observed experimentally by 
Calvetti et al. (HHI on a Schneebeli assembly |10j prepared using something like the 
GG procedure, and also by Radjai et al. [JTJ H21 HH] with discrete element method 
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simulations (Contact Dynamics). In the same vein, Geng et al. [43] found that 
a raining procedure for building 2D heaps of photoelastic particles resulted in more 
nearly axisymmetric distributions of contacts than did a point-source pouring method. 
Interestingly, in these latter experiments, the underlying symmetry was hexagonal rather 
than square. At the system size level, the difference is less impressive: we have computed 
the solid fraction <p and found 4>rl — 0.820(2) and 4>gg — 0.812(4) respectively. The 
average numbers of contacts per grain, Z, are also rather close: Zrl — 3.39(1) and 
Z GG = 3.517(8). 

3.2. The response function 

After the deposition procedure converges to a static equilibrium, which typically takes 
10 6 MD time steps, an additional force F localized on a single grain is added at the top 
of the layer. Its value is typically of the order of a few times the weight of a grain. As 
we emphasize below, F$ can be either vertical or inclined at a 45°. This force is applied 
progressively over 10 4 MD time steps, and we checked that this overloading does not 
lead to important rearrangements of the packing. More precisely, no contacts are either 
broken or created, and the few observed sliding events do not push the response out of 
a linear and reversible regime. Figure shows a map of the contact forces in response 
to this extra force - the force network before the overloading was subtracted. In order 
to get the averaged stress response function, several tens of such independent overloads 
were performed on assemblies of 3600 particles. The stress components a zz and o~ xz are 
measured at the bottom of the system in the following manner. Like the experimental 
pressure probe of size Z, the stress component a az at position x is obtained as the sum 
over the a-components of the forces exerted on the bottom grains located at horizontal 
distances between x — 1/2 and x + £/2 from the overload point. This sum is divided by 
the coarse graining length, £, and by the number of samples M. It is also normalized by 
the force F . Finally: 

, , 1 



cr nz (x 



E mi (i) 



We checked that this procedure yields results consistent with a more general definition 
of the stress tensor a a /3 [El 1221 EE GUI EZ| • We note that the stress components become 
independent of the coarse graining length only when i is large enough compared to the 
mean grain diameter d Here we used i ~ 7.5d, which is about the crossover value. 

Figures |H1 and El show the different stress response profiles. They exhibit a single 
broad peak. Although not shown here, and in accordance with the experimental data 
presented in the previous sections, these profiles also scale linearly with the layer 
thickness. The o zz profiles measured from the RL or the GG procedures are distinct but 
close. Interestingly, the difference between the two preparations is much more visible 
on the a xz plots. In particular, the amplitude of the shear maximum is typically 10% 
larger in the more anisotropic case (GG). 
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Figure 9. Shear stress a xz profiles in response to a vertical (left) or a 45° inclined 
(right) top force. The difference between the two preparation procedures is much more 
visible than in the data for a zz . In particular, the amplitude of the shear maximum is 
~ 10% larger for the more anisotropic case. 



3.3. Comparison with isotropic elasticity 

As in the loose and dense cases of the previous section, both rain-like and grain-by- 
grain preparations have the vertical axis as an axis of symmetry. It is not clear a 
priori that any piling obtained under such conditions would ever be isotropic. From a 
microscopic point of view, an inspection of the contact distribution (figure EJ) suggests 
that that the RL procedure is more likely to yield an isotropic effective medium than 
the GG deposition procedure. The question is how much are the texture properties 
reflected in macroscopic features. Here, we simply try to adjust the stress response 
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Figure 10. Comparison of the stress profiles with isotropic elasticity. The fit for both 
cases of a vertical and an inclined overload Fq is quite good. These plots are for the 
rain-like preparation of the granular layer. 
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Figure 11. Equivalent of figure ITU1 for the case of the grain- by-grain preparation. 
Here, elastic predictions of isotropic elasticity for both vertical and non-normal extra 
forces are much less good. 



functions predicted by a isotropic elasticity, and we seek to identify which features of 
these functions would be most sensitive for identifying a departure from isotropy. 

In the same vein as the experiments discussed above, j20j, we performed fits of 
the numerical profiles by minimizing the quadratic distance between the data and the 
theoretical predictions, including the weight of the uncertainties. Recall that the fitting 
parameter is the Poisson coefficient v. As shown in figure EB the RL profiles can be well 
reproduced with v = 0.28(2) (best fit) when the extra force F is vertical. This value 
of the Poisson ratio also fits well the stresses in response to an overload inclined at 45°. 
In the case of the GG preparation ( figure ITTj) . the best fit of the stress responses to a 
vertical F {y = 0.25(2)) is not as well adjusted to the numerical data. The discrepancy 
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between isotropic elasticity and the GG data is also clearly visible for the inclined F 
situation. Using a range of directions for the applied extra force provides a good test 
to identify a medium showing weak anisotropy such as occurs for preparations under 
gravity. 

4. Discussion and conclusion 

Here we presented a summary of recent results based on the measurement of stress 
responses to localized forces. These results are basically consistent with an elastic 
description. We have shown experimentally for 2D model granular assemblies and for 
3D sand piles that the response function is a useful probe of constitutive properties of 
the granular material. Indeed, this type of probe is a non intrusive method that could 
replace or supplement standard acoustic techniques, since it is well adapted to weakly 
cohesive material under low confining pressures where sound propagation methods are 
more difficult to implement. 

We have shown that the response function is sensitive to small differences in the 
preparation history and also that its shape may change significantly after a quasi-static 
deformation. Using numerical simulations, we associate these macroscopic differences 
to different microscopic features which appear in the angular distribution of contacts. 
Nevertheless, exactly how these features would quantitatively emerge at the macroscopic 
level and how they are likely to change due to an external loading path is still a 
challenging and unanswered question. 

Finally, we propose a new path of systematic measurements where the normal and 
shear stress profiles are obtained for various directions of the applied localized force. 
A collection of such data would yield crucial information on the structure evolution 
for a given loading history and, therefore, provide a severe test of any model such as 
anisotropic elasticity models. 
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